{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# System ready for equilibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from htmd.ui import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "path=mkdtemp() + '/equilibrate'\n",
    "copytree(home()+'/data/equilibrate', path )\n",
    "chdir(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These steps show how to prepare your system for equilibration. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check for the acemd protocols available within your HTMD distribution,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['equilibration.py']\n"
     ]
    }
   ],
   "source": [
    "Acemd.protocols()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "set numsteps 2500000\n",
      "set temperature 300\n",
      "proc calcforces_init {} {\n",
      "  berendsenpressure  off}\n",
      "proc calcforces {} {\t\n",
      "  global numsteps\n",
      "  set step [ getstep ]\n",
      "  if { $step > 500 } {\n",
      "    berendsenpressure  on\n",
      "  } else {\n",
      "    berendsenpressure  off}\n",
      "  if { $step > [expr $numsteps/2] } {\n",
      "    constraintscaling 0\n",
      "  } else {\n",
      "    constraintscaling [expr 1 + $step*(0.1-1)*2/$numsteps]}\n",
      "}\n",
      "proc calcforces_endstep { } { }\n",
      "1-4scaling\t1.0\n",
      "cutoff\t9\n",
      "hydrogenscale\t4\n",
      "consref\tstructure.pdb\n",
      "tclforces\ton\n",
      "constraints\ton\n",
      "langevindamping\t1\n",
      "switching\ton\n",
      "minimize\t500\n",
      "structure\tstructure.psf\n",
      "restartfreq\t5000\n",
      "temperature\t$temperature\n",
      "exclude\tscaled1-4\n",
      "berendsenpressuretarget\t1.01325\n",
      "constraintscaling\t1.0\n",
      "pme\ton\n",
      "switchdist\t7.5\n",
      "xtcfile\toutput.xtc\n",
      "langevintemp\t$temperature\n",
      "rigidbonds\tall\n",
      "fullelectfrequency\t2\n",
      "berendsenpressure\ton\n",
      "outputname\toutput\n",
      "langevin\ton\n",
      "restart\ton\n",
      "xtcfreq\t25000\n",
      "energyfreq\t1000\n",
      "timestep\t4\n",
      "berendsenpressurerelaxationtime\t800\n",
      "run\t$numsteps\n",
      "pmegridspacing\t1.0\n",
      "parameters\tparameters\n",
      "coordinates\tstructure.pdb\n",
      "\n"
     ]
    }
   ],
   "source": [
    "acemd = Acemd('equilibration.py')\n",
    "acemd.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For use of constraints during the equilibration set occupancy and beta columns of protein and ligand to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s = Molecule('structure.pdb')\n",
    "s.set('occupancy',0)\n",
    "s.set('beta',0)\n",
    "s.set('beta',1,sel='segid L and noh')\n",
    "s.set('beta',1,sel='segid P and noh')\n",
    "s.write('structure.pdb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate ans set the size of the system periodic box,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "box = amax(s.get('coords','water'),axis=0)-amin(s.get('coords','water'),axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "acemd.celldimension = str(box[0])+' '+str(box[1])+' '+str(box[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Provide here the name of the directory you used to build your system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "acemd.load('./')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the input folder ready for running an equilibration of your system using ACEMD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "acemd.save('./equil/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start the equilibration of this system, just run ACEMD from the ./equil directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "loc = LocalGPUQueue()\n",
    "loc.submit('./equil/')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}